############################
# Replication code for:
# Awareness of Executive Interference and the Demand for Judicial Independence: Evidence from Four Constitutional Courts
############################

### This code replicates the following results:

# Main text
  # Figures 1 and 2 (V-Dem indicators): Lines 878-1062
  # Figure 3 (Distribution of Awareness): Lines 71-95
  # Figure 4 (Distribution of Perceived Executive Influence): Lines 97-121
  # Figure 5 (Distribution of Demand for Judicial Independence): Lines 123-150
  # Table 2 (Determinants of Executive Influence): Lines 159-219
  # Table 3 (Determinants of Demand for Judicial Independence): Lines 274-343
  # Figure 6: Lines 345-489

# Appendix
  # Tables B2-B5 (Variable Statistics): Lines 697-767
  # Table C1 (Full regression results): Lines 492-523
  # Table C2 (Bivariate models): Lines 525-590
  # Table C3 (Ordered logit models): Lines 592-695
  # Table C4 (Pooled models): Lines 222-271
  # Table D1 (Construction of Diffuse Support Scale): Lines 769-875


# Required packages
library(sandwich)
library(lmtest)
library(marginaleffects)
library(MASS)
library(stargazer)
library(xtable)
library(psych)
library(performance)

# Set directory
setwd()

# Clean environment
rm(list=ls())

load("data/data_JLC.RData")

names(df)

# Country samples
df_US <- df[df$country == "U.S." & complete.cases(df), ]
df_DE <- df[df$country == "Germany" & complete.cases(df), ]
df_HU <- df[df$country == "Hungary" & complete.cases(df), ]
df_PL <- df[df$country == "Poland" & complete.cases(df), ]

# Pooled sample
df_all <- df[complete.cases(df), ]

##################
# Distribution of main variables
# Main text - Figures 3-5
##################

## Grouped bar charts
colors <- viridis::viridis(4)
col3 <- c(colors[2], adjustcolor(colors[2], alpha.f = 0.66),
          adjustcolor(colors[2], alpha.f = 0.33))
col4 <- c(colors[2], adjustcolor(colors[2], alpha.f = 0.75),
          adjustcolor(colors[2], alpha.f = 0.46875),
          adjustcolor(colors[2], alpha.f = 0.1875))
# Check colors
barplot(c(1:3), col=col3)
barplot(c(1:4), col=col4)

### Awareness (Figure 3)
table(df_US$aware_ct, useNA = "always")

aware <- matrix(NA, ncol = 4, nrow=4)
colnames(aware) <- c("U.S.", "Germany", "Hungary", "Poland")
rownames(aware) <- c("Very aware","Somewhat aware","Not very aware", "Never heard of it")

aware[,1] <- as.numeric(table(df_US$aware_ct)/nrow(df_US)*100)
aware[,2] <- as.numeric(table(df_DE$aware_ct)/nrow(df_DE)*100)
aware[,3] <- as.numeric(table(df_HU$aware_ct)/nrow(df_HU)*100)
aware[,4] <- as.numeric(table(df_PL$aware_ct)/nrow(df_PL)*100)
aware

pdf(file = "fig/JLC-fig3.pdf", width = 8, height = 4)
# .tiff file
#tiff(file = "fig/awareness-barplot.tif", width = 8, height = 4, units = "in", res = 300)
par(mar=c(3, 3.35, 1, 10), xpd=TRUE)
barplot(aware, col=rev(col4), xlab="", ylab = "", 
        main = "", cex.names = 1.1)

mtext(side = 2, "Percentage of respondents", line = 2.35)

legend("topright", legend = rownames(aware), inset=c(-0.39,0),
       fill = col4, bty = "n", y.intersp = 1, cex = 1.3)
dev.off()

### Perceived Executive Influence (Figure 4)
table(df_US$exec_infl, useNA = "always")

perc <- matrix(NA, ncol = 4, nrow=4)
colnames(perc) <- c("U.S.", "Germany","Hungary","Poland")
rownames(perc) <- c("Great extent","Moderate extent","Small extent", "Not at all")

perc[,1] <- as.numeric(table(df_US$exec_infl)/nrow(df_US)*100)
perc[,2] <- as.numeric(table(df_DE$exec_infl)/nrow(df_DE)*100)
perc[,3] <- as.numeric(table(df_HU$exec_infl)/nrow(df_HU)*100)
perc[,4] <- as.numeric(table(df_PL$exec_infl)/nrow(df_PL)*100)
perc

pdf(file = "fig/JLC-fig4.pdf", width = 8, height = 4)
# .tiff file
#tiff(file = "fig/exec-infl-barplot.tif", width = 8, height = 4, units = "in", res = 300)
par(mar=c(3, 3.35, 1, 10), xpd=TRUE)
barplot(perc, col=rev(col4), xlab="", ylab = "",
        main = "", cex.names = 1.1)

mtext(side = 2, "Percentage of respondents", line = 2.35)

legend("topright", legend = rownames(perc), inset=c(-0.38,0),
       fill = col4, bty = "n", y.intersp = 1, cex = 1.3)
dev.off()

## Demand for Judicial Independence (Figure 5)
table(df_US$demand_ord, useNA = "always")
table(df_DE$demand_ord, useNA = "always")
table(df_HU$demand_ord, useNA = "always")
table(df_PL$demand_ord, useNA = "always")

demand <- matrix(NA, ncol = 4, nrow=3)
colnames(demand) <- c("U.S.", "Germany","Hungary","Poland")
rownames(demand) <- c("Too much", "About right", "Too little")

demand[,1] <- as.numeric(table(df_US$demand_ord)/nrow(df_US)*100)
demand[,2] <- as.numeric(table(df_DE$demand_ord)/nrow(df_DE)*100)
demand[,3] <- as.numeric(table(df_HU$demand_ord)/nrow(df_HU)*100)
demand[,4] <- as.numeric(table(df_PL$demand_ord)/nrow(df_PL)*100)
demand

pdf(file = "fig/JLC-fig5.pdf", width = 8, height = 4)
# .tiff file
#tiff(file = "fig/demand-barplot.tif", width = 8, height = 4, units = "in", res = 300)
par(mar=c(3, 3.35, 1, 10), xpd=TRUE)
barplot(demand, col=rev(col3), xlab = "", ylab = "",
        main = "", cex.names = 1.1)

mtext(side = 2, "Percentage of respondents", line = 2.35)

legend("topright", legend = rownames(demand), inset=c(-0.3,0),
       fill = col3, bty = "n", y.intersp = 1, cex = 1.3)
dev.off()


###############################
# Determinants of perceived executive influence
# Main text - Table 2
# Appendix - Table C4
###############################

### Main text Table 2

# Bivariate models
aware_biv <- as.formula("exec_infl ~ aware_ct")

ols_US_biv <- lm(aware_biv, data = df_US)
se_US_biv <- coeftest(ols_US_biv, vcov = vcovHC(ols_US_biv, type = "HC3"))[,2]

ols_DE_biv <- lm(aware_biv, data = df_DE)
se_DE_biv <- coeftest(ols_DE_biv, vcov = vcovHC(ols_DE_biv, type = "HC3"))[,2]

ols_HU_biv <- lm(aware_biv, data = df_HU)
se_HU_biv <- coeftest(ols_HU_biv, vcov = vcovHC(ols_HU_biv, type = "HC3"))[,2]

ols_PL_biv <- lm(aware_biv, data = df_PL)
se_PL_biv <- coeftest(ols_PL_biv, vcov = vcovHC(ols_PL_biv, type = "HC3"))[,2]

# Models with controls
aware_controls <- as.formula("exec_infl ~ aware_ct + df_support + perf_ct + copartisan + ideo10 + pol_interest + strong_lead + know_ct + age + gender")

ols_US_controls <- lm(aware_controls, data = df_US)
se_US_controls <- coeftest(ols_US_controls, vcov = vcovHC(ols_US_controls, type = "HC3"))[,2]

ols_DE_controls <- lm(aware_controls, data = df_DE)
se_DE_controls <- coeftest(ols_DE_controls, vcov = vcovHC(ols_DE_controls, type = "HC3"))[,2]

ols_HU_controls <- lm(aware_controls, data = df_HU)
se_HU_controls <- coeftest(ols_HU_controls, vcov = vcovHC(ols_HU_controls, type = "HC3"))[,2]

ols_PL_controls <- lm(aware_controls, data = df_PL)
se_PL_controls <- coeftest(ols_PL_controls, vcov = vcovHC(ols_PL_controls, type = "HC3"))[,2]

# Prepare Table 2
labels <- c("Awareness", "Diffuse Supp.", "Specific Support", "Gov. Supporter",
            "Ideology", "Pol. Interest", "Strong Leader", "Knowledge Court", "Age", "Gender")

# Table 2
tab2 <- stargazer(ols_US_biv, ols_US_controls, 
        ols_DE_biv, ols_DE_controls,
        ols_HU_biv, ols_HU_controls, 
        ols_PL_biv, ols_PL_controls, 
        se = list(se_US_biv, se_US_controls, 
          se_DE_biv, se_DE_controls, 
          se_HU_biv, se_HU_controls, 
          se_PL_biv, se_PL_controls),
        keep.stat = c("n","adj.rsq"),
        star.char = c("+", "*", "**", "***"),
        star.cutoffs = c(.1, .05, .01, .001), 
        no.space=TRUE, align=TRUE, table.placement = "H", 
        title="Determinants of Perceived Executive Influence", 
        label = "tab:exec-infl",
        dep.var.labels.include = FALSE, header = FALSE,
        dep.var.caption = "",
        column.labels = c("\\textbf{US}", "\\textbf{DE}",
                        "\\textbf{HU}", "\\textbf{PL}"),
        column.separate = c(2,2,2,2),
        covariate.labels = labels,
        column.sep.width = "-30pt",
        notes.append=FALSE, notes.align = "l", notes.label = "",
        notes = c("\\footnotesize \\parbox[t]{\\textwidth}{\\textit{Note:} Robust (heteroskedasticity-consistent) standard errors in parentheses. \\hfill $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}"))
cat(tab2, sep = '\n', file = "tab/JLC-tab2.tex")

### Pooled models - Appendix Table C4
  # Baseline: US
f_pooled <- as.formula("exec_infl ~ aware_ct*country_base + df_support + perf_ct + copartisan + ideo10 + pol_interest + strong_lead + know_ct + age + gender")
df_all$country_base <- relevel(df_all$country, ref = "U.S.")
pooled_exec_US <- lm(f_pooled, data = df_all)
  # Baseline: DE
df_all$country_base <- relevel(df_all$country, ref = "Germany")
pooled_exec_DE <- lm(f_pooled, data = df_all)
  # Baseline: HU
df_all$country_base <- relevel(df_all$country, ref = "Hungary")
pooled_exec_HU <- lm(f_pooled, data = df_all)
  # Baseline: PL
df_all$country_base <- relevel(df_all$country, ref = "Poland")
pooled_exec_PL <- lm(f_pooled, data = df_all)

se_pooled_exec_US <- coeftest(pooled_exec_US, vcovCL(pooled_exec_US, cluster = ~ country))[,2]
se_pooled_exec_DE <- coeftest(pooled_exec_DE, vcovCL(pooled_exec_DE, cluster = ~ country))[,2]
se_pooled_exec_HU <- coeftest(pooled_exec_HU, vcovCL(pooled_exec_HU, cluster = ~ country))[,2]
se_pooled_exec_PL <- coeftest(pooled_exec_PL, vcovCL(pooled_exec_PL, cluster = ~ country))[,2]

labels <- c("Awareness", 
            "Germany (DE)", "Hungary (HU)", "Poland (PL)", "United States (US)",
            "Awareness $\\times$ DE", "Awareness $\\times$ HU", 
            "Awareness $\\times$ PL", "Awareness $\\times$ US")

# Appendix Table C4
tabC4 <- stargazer(pooled_exec_US, pooled_exec_DE, 
          pooled_exec_HU, pooled_exec_PL,
          se = list(se_pooled_exec_US, se_pooled_exec_DE, 
                    se_pooled_exec_HU, se_pooled_exec_PL),
          keep.stat = c("n","adj.rsq"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001), 
          no.space=TRUE, align=TRUE, table.placement = "H", 
          title="Determinants of Executive Influence (Pooled Models)", 
          label = "tab:pooled-awareness", dep.var.caption = "",
          dep.var.labels.include = FALSE, header = FALSE,
          column.labels = c("\\textbf{US}", "\\textbf{DE}",
                            "\\textbf{HU}", "\\textbf{PL}"),
          keep = c("aware_ct", 
                   "country_baseGermany", 
                   "country_baseHungary", 
                   "country_basePoland", 
                   "country_baseU.S.",
                   "Constant"),
          covariate.labels = labels,
          column.sep.width = "-10pt",
          notes.append=FALSE, notes.align = "l", notes.label = "",
          notes = c("\\footnotesize \\parbox[t]{\\textwidth}{\\textit{Note:} Robust standard errors clustered by country in parentheses. \\hfill $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}"))
cat(tabC4, sep = '\n', file = "tab/JLC-tabC4.tex")


###############################
# Determinants of demand for judicial independence
# Main text - Table 3
# Main text - Figure 6
###############################

### Main text Table 3

# Additive models
demand_add <- as.formula("demand ~ aware_ct + exec_infl + df_support + perf_ct + copartisan + ideo10 + pol_interest + strong_lead + know_ct + age + gender")

ols_US <- lm(demand_add, data = df_US)
vcovUS <- vcovHC(ols_US, type = "HC3")
se_US <- coeftest(ols_US, vcov = vcovHC(ols_US, type = "HC3"))[,2]

ols_DE <- lm(demand_add, data = df_DE)
vcovDE <- vcovHC(ols_DE, type = "HC3")
se_DE <- coeftest(ols_DE, vcov = vcovHC(ols_DE, type = "HC3"))[,2]

ols_HU <- lm(demand_add, data = df_HU)
vcovHU <- vcovHC(ols_HU, type = "HC3")
se_HU <- coeftest(ols_HU, vcov = vcovHC(ols_HU, type = "HC3"))[,2]

ols_PL <- lm(demand_add, data = df_PL)
vcovPL <- vcovHC(ols_PL, type = "HC3")
se_PL <- coeftest(ols_PL, vcov = vcovHC(ols_PL, type = "HC3"))[,2]

# Interaction models
demand_int <- as.formula("demand ~ aware_ct*exec_infl + df_support + perf_ct + copartisan + ideo10 + pol_interest + strong_lead + know_ct + age + gender")

ols_US1 <- lm(demand_int, data = df_US)
vcovUS1 <- vcovHC(ols_US1, type = "HC3")
se_US1 <- coeftest(ols_US1, vcov = vcovHC(ols_US1, type = "HC3"))[,2]

ols_DE1 <- lm(demand_int, data = df_DE)
vcovDE1 <- vcovHC(ols_DE1, type = "HC3")
se_DE1 <- coeftest(ols_DE1, vcov = vcovHC(ols_DE1, type = "HC3"))[,2]

ols_HU1 <- lm(demand_int, data = df_HU)
vcovHU1 <- vcovHC(ols_HU1, type = "HC3")
se_HU1 <- coeftest(ols_HU1, vcov = vcovHC(ols_HU1, type = "HC3"))[,2]

ols_PL1 <- lm(demand_int, data = df_PL)
vcovPL1 <- vcovHC(ols_PL1, type = "HC3")
se_PL1 <- coeftest(ols_PL1, vcov = vcovHC(ols_PL1, type = "HC3"))[,2]

labels <- c("Awareness", "Perceived Exec. Infl.", "Awareness $\\times$ Exec. Infl.")

# Table 3 (short table only displaying main variables)
tab3 <- stargazer(ols_US, ols_DE, ols_HU, ols_PL,
          ols_US1, ols_DE1, ols_HU1, ols_PL1,
          se = list(se_US, se_DE, se_HU, se_PL,
                    se_US1, se_DE1, se_HU1, se_PL1),
          keep.stat = c("n","adj.rsq"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001), 
          no.space=TRUE, align=TRUE, table.placement = "H", 
          title="Determinants of Demand for Judicial Independence", 
          label = "tab:3",
          dep.var.labels.include = FALSE, header = FALSE,
          dep.var.caption = "",
          column.labels = c("\\textbf{Additive Model}",
                            "\\textbf{Interaction Model}"),
          column.separate = c(4, 4),
          keep = c("aware_ct", "exec_infl", "Constant"),
          covariate.labels = labels,
          column.sep.width = "-32pt",
          notes.append=FALSE, notes.align = "l", notes.label = "",
          notes = c("\\footnotesize \\parbox[t]{\\textwidth}{\\textit{Note:} Controls:	\\textit{Diffuse Support}, \\textit{Specific Support}, \\textit{Partisanship (Gov. Supporter)}, \\textit{Ideology}, \\textit{Political Interest}, \\textit{Strong Leader}, \\textit{Knowledge of the Court}, \\textit{Age}, and \\textit{Gender}. See Appendix C for full regression results. Robust (heteroskedasticity-consistent) standard errors in parentheses. \\hfill $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}"))
cat(tab3, sep = '\n', file = "tab/JLC-tab3.tex")

### Main text - Figure 6

# Compute coefficient estimates and SEs for figure 6 (left panel)
estimates <- rbind(ols_US$coefficients["exec_infl"], ols_DE$coefficients["exec_infl"], 
                   ols_HU$coefficients["exec_infl"], ols_PL$coefficients["exec_infl"])
CI <- rbind(
  estimates[1] + (se_US[3] * qt(c(0.035,0.975),ols_US$df)),
  estimates[2] + (se_DE[3] * qt(c(0.035,0.975),ols_DE$df)), 
  estimates[3] + (se_HU[3] * qt(c(0.035,0.975),ols_HU$df)), 
  estimates[4] + (se_PL[3] * qt(c(0.035,0.975),ols_PL$df))
)

plot_coef <- cbind(estimates, CI)
colnames(plot_coef) <- c("estimate", "lower", "upper")
rownames(plot_coef) <- c("US","DE", "HU", "PL")
plot_coef <- as.data.frame(plot_coef)
plot_coef$country <- c("US","DE", "HU", "PL")

# Compute marginal effects for figure 6 (right panel)
meff_US1 <- cbind(data.frame(avg_comparisons(ols_US1, variable = "aware_ct", by = "exec_infl", 
                                             vcov = vcovUS1))[c(1:7,9:10)], country = "US")
meff_DE1 <- cbind(data.frame(avg_comparisons(ols_DE1, variable = "aware_ct", by = "exec_infl", 
                                             vcov = vcovDE1))[c(1:7,9:10)], country = "DE")
meff_HU1 <- cbind(data.frame(avg_comparisons(ols_HU1, variable = "aware_ct", by = "exec_infl", 
                                             vcov = vcovHU1))[c(1:7,9:10)], country = "HU")
meff_PL1 <- cbind(data.frame(avg_comparisons(ols_PL1, variable = "aware_ct", by = "exec_infl", 
                                             vcov = vcovPL1))[c(1:7,9:10)], country = "PL")

plot_meff1 <- rbind(meff_US1, meff_DE1, meff_HU1, meff_PL1)


pdf(file = "fig/JLC-fig6.pdf", width = 9, height = 4.25) # height of plot (inches)
# .tiff file
#tiff(file = "fig/JLC-fig.tif", width = 9, height = 4.25, units = "in", res = 300)
par(mar=c(3, 2.35, 0.5, 1.5), # c(bottom, left, top, right)
    oma=c(0.5, 0.5, 0.5, 0.5))
layout(matrix(1:2, ncol=2), # set up panel
       widths=c(0.345, 0.655)) # relative widths

# Left panel

# Set panel features
ylim <- c(-0.05, 1.05)
yticks <- c(0, 0.33, 0.66, 1) # location of y labels
ylabels <- plot_coef$country

plot(1, 1, type = "n", yaxt = "n", xaxt = "n",
     xlim = c(0, max(plot_coef$lower, plot_coef$upper)), 
     ylim = rev(ylim),
     xlab = "", 
     ylab = "", main = "")

axis(side = 2, at=yticks, las=1, labels = ylabels, tick = TRUE, 
     cex.axis = 1.3)

axis(side = 1, at=seq(0,0.3,0.1), las=1, labels = FALSE, cex.axis = 1.1)

text(x = seq(0,0.3,0.1), # horizontal position of x labels
     y = 1.16,
     labels = seq(0,0.3,0.1),
     adj = 0.5, # alignment (0=left, 0.5=center, 1=right)
     xpd = NA,
     cex = 1)

mtext("Coefficient on Perceived Exec. Influence", side = 1, line = 2.35, 
      cex = 1.1)

segments(x0=plot_coef$lower, x1=plot_coef$upper, 
         y0=yticks, y1=yticks, lwd = 3, col = "black")

points(x = plot_coef$estimate, y = yticks, pch = 19, cex = 1.5, col = "black")

abline(v = 0, col = "red", lty = 2, lwd = 1.5)

# Right panel

# Set panel features
xlim <- c(-0.03, 1.03)
xticks <- c(0, 0.33, 0.66, 1) # location of x labels
# location of point estimates (separation between country estimates)
pticks <- cbind(c(-0.03, -0.01, 0.01, 0.03),
                c(0.30, 0.32, 0.34, 0.36),
                c(0.63, 0.65, 0.67, 0.69),
                c(0.97, 0.99, 1.01, 1.03))
xlabels <- c("Not \nat all", "Small\nextent", "Moderate\nextent", "Great\nextent")
legend <- unique(plot_meff1$country)

par(mar=c(4, 3, 0.5, 0.5))

plot(1, 1, type = "n", xaxt = "n",
     ylim = c(min(plot_meff1$conf.low, plot_meff1$conf.high)-0.01, 
              max(plot_meff1$conf.low, plot_meff1$conf.high)+0.01), 
     xlim = xlim,
     ylab = "", main = "",
     xlab = "")
# y label closer to axis
mtext("Marginal Effect of Awareness", side = 2, line = 2.35)

mtext("Perceived Executive Influence", side = 1, line = 3.25, cex = 1.1)

# x axis
axis(side = 1, at=xticks, las=1, labels = xlabels, tick = FALSE)

# US
segments(y0=plot_meff1$conf.low[plot_meff1$country=="US"], 
         y1=plot_meff1$conf.high[plot_meff1$country=="US"], 
         x0=pticks[1,], x1=pticks[1,], lwd = 2, col = "black")
lines(y = plot_meff1$estimate[plot_meff1$country=="US"], x = pticks[1,], lwd = 2.5, lty = 2,
      bg = "black")
points(y = plot_meff1$estimate[plot_meff1$country=="US"], x = pticks[1,], pch = 21, cex = 1.3, 
       bg = "black", col = "black")

# DE
segments(y0=plot_meff1$conf.low[plot_meff1$country=="DE"], 
         y1=plot_meff1$conf.high[plot_meff1$country=="DE"], 
         x0=pticks[2,], x1=pticks[2,], lwd = 2, col = "black")
lines(y = plot_meff1$estimate[plot_meff1$country=="DE"], x = pticks[2,], lwd = 2.5, lty = 32,
      bg = "black")
points(y = plot_meff1$estimate[plot_meff1$country=="DE"], x = pticks[2,], pch = 24, cex = 1.3, 
       bg = "black", col = "black")

# HU
segments(y0=plot_meff1$conf.low[plot_meff1$country=="HU"], 
         y1=plot_meff1$conf.high[plot_meff1$country=="HU"], 
         x0=pticks[3,], x1=pticks[3,], lwd = 2, col = "grey")
lines(y = plot_meff1$estimate[plot_meff1$country=="HU"], x = pticks[3,], lwd = 2.5, lty = 1,
      col = "grey")
points(y = plot_meff1$estimate[plot_meff1$country=="HU"], x = pticks[3,], pch = 21, cex = 1.3, 
       bg = "grey", col = "grey")
# PL
segments(y0=plot_meff1$conf.low[plot_meff1$country=="PL"], 
         y1=plot_meff1$conf.high[plot_meff1$country=="PL"], 
         x0=pticks[4,], x1=pticks[4,], lwd = 2, col = "grey")
lines(y = plot_meff1$estimate[plot_meff1$country=="PL"], x = pticks[4,], lwd = 2.5, lty = 1,
      col = "grey")
points(y = plot_meff1$estimate[plot_meff1$country=="PL"], x = pticks[4,], pch = 24, cex = 1.3, 
       bg = "grey", col = "grey")

legend("topleft", legend, lwd = 1.5, lty = c(2,2,1,1), pch = c(21,24,21,24), 
       col = c("black", "black", "grey", "grey"), ncol = 2,
       pt.bg = c("black", "black", "grey", "grey"))

abline(h = 0, col = "red", lty = 2, lwd = 2)

dev.off()


###############################
# Determinants of demand for judicial independence
# Appendix - Tables C1-C3
###############################

### Appendix Table C1 (full table)

labels_full <- c("Awareness", "Perceived Exec. Infl.", "Awareness $\\times$ Exec. Infl.", 
                 "Diffuse Supp.", "Specific Support", "Gov. Supporter", "Ideology", 
                 "Pol. Interest", "Strong Leader", "Knowledge Court", "Age", "Gender")

tabC1 <- stargazer(ols_US, ols_DE, ols_HU, ols_PL,
          ols_US1, ols_DE1, ols_HU1, ols_PL1,
          se = list(se_US, se_DE, se_HU, se_PL,
               se_US1, se_DE1, se_HU1, se_PL1),
          keep.stat = c("n","adj.rsq"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001), 
          no.space=TRUE, align=TRUE, table.placement = "H", 
          title="Determinants of Demand for Judicial Independence (OLS Models)", 
          label = "tab:C1",
          dep.var.labels.include = FALSE, header = FALSE,
          dep.var.caption = "",
          column.labels = c("\\textbf{Additive Model}",
                            "\\textbf{Interaction Model}"),
          column.separate = c(4, 4),
          order = c(1:2, 12, 2:11,13),
          covariate.labels = labels_full,
          column.sep.width = "-32pt",
          notes.append=FALSE, notes.align = "l", notes.label = "",
          notes = c("\\footnotesize \\parbox[t]{\\textwidth}{\\textit{Note:} Robust (heteroskedasticity-consistent) standard errors in parentheses. \\hfill $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}"))
cat(tabC1, sep = '\n', file = "tab/JLC-tabC1.tex")

### Appendix Table C2 (bivariate models)

# Bivariate models - Awareness
demand_aware <- as.formula("demand ~ aware_ct")

ols_US_aware <- lm(demand_aware, data = df_US)
vcovUS_aware <- vcovHC(ols_US_aware, type = "HC3")
se_US_aware <- coeftest(ols_US_aware, vcov = vcovHC(ols_US_aware, type = "HC3"))[,2]

ols_DE_aware <- lm(demand_aware, data = df_DE)
vcovDE_aware <- vcovHC(ols_DE_aware, type = "HC3")
se_DE_aware <- coeftest(ols_DE_aware, vcov = vcovHC(ols_DE_aware, type = "HC3"))[,2]

ols_HU_aware <- lm(demand_aware, data = df_HU)
vcovHU_aware <- vcovHC(ols_HU_aware, type = "HC3")
se_HU_aware <- coeftest(ols_HU_aware, vcov = vcovHC(ols_HU_aware, type = "HC3"))[,2]

ols_PL_aware <- lm(demand_aware, data = df_PL)
vcovPL_aware <- vcovHC(ols_PL_aware, type = "HC3")
se_PL_aware <- coeftest(ols_PL_aware, vcov = vcovHC(ols_PL_aware, type = "HC3"))[,2]

# Bivariate models - Executive Influence
demand_exec <- as.formula("demand ~ exec_infl")

ols_US_exec <- lm(demand_exec, data = df_US)
vcovUS_exec <- vcovHC(ols_US_exec, type = "HC3")
se_US_exec <- coeftest(ols_US_exec, vcov = vcovHC(ols_US_exec, type = "HC3"))[,2]

ols_DE_exec <- lm(demand_exec, data = df_DE)
vcovDE_exec <- vcovHC(ols_DE_exec, type = "HC3")
se_DE_exec <- coeftest(ols_DE_exec, vcov = vcovHC(ols_DE_exec, type = "HC3"))[,2]

ols_HU_exec <- lm(demand_exec, data = df_HU)
vcovHU_exec <- vcovHC(ols_HU_exec, type = "HC3")
se_HU_exec <- coeftest(ols_HU_exec, vcov = vcovHC(ols_HU_exec, type = "HC3"))[,2]

ols_PL_exec <- lm(demand_exec, data = df_PL)
vcovPL_exec <- vcovHC(ols_PL_exec, type = "HC3")
se_PL_exec <- coeftest(ols_PL_exec, vcov = vcovHC(ols_PL_exec, type = "HC3"))[,2]

labels_biv <- c("Awareness", "Perceived Exec. Infl.")

tabC2 <- stargazer(ols_US_aware,ols_US_exec, 
          ols_DE_aware, ols_DE_exec,
          ols_HU_aware, ols_HU_exec,
          ols_PL_aware, ols_PL_exec,
          se = list(se_US_aware, se_US_exec,
                    se_DE_aware, se_DE_exec,
                    se_HU_aware, se_HU_exec,
                    se_PL_aware,se_PL_exec),
          keep.stat = c("n","adj.rsq"),
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001), 
          no.space=TRUE, align=TRUE, table.placement = "H", 
          title="Determinants of Demand for Judicial Independence (Bivariate Models)", 
          label = "tab:ols-demand-bivariate",
          dep.var.labels.include = FALSE, header = FALSE,
          dep.var.caption = "",
          column.labels = c("\\textbf{US}", "\\textbf{DE}", 
                            "\\textbf{HU}", "\\textbf{PL}"),
          column.separate = c(2,2,2,2),
          covariate.labels = labels_biv,
          column.sep.width = "-35pt",
          notes.append=FALSE, notes.align = "l", notes.label = "",
          notes = c("\\footnotesize \\parbox[t]{\\textwidth}{\\textit{Note:} Robust (heteroskedasticity-consistent) standard errors in parentheses. \\hfill $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}"))
cat(tabC2, sep = '\n', file = "tab/JLC-tabC2.tex")

### Appendix Table C3 (ordered logit models)

table(df_US$demand_ord)
levels(df_US$demand_ord)

# Additive model
demand.ord <- as.formula("demand_ord ~ aware_ct + exec_infl + df_support + perf_ct + copartisan + ideo10 + pol_interest + strong_lead + know_ct + age + gender")

olog_US <- polr(demand.ord, Hess = TRUE, method = "logistic", data=df_US)

olog_DE <- polr(demand.ord, Hess = TRUE, method = "logistic", data=df_DE)

olog_HU <- polr(demand.ord, Hess = TRUE, method = "logistic", data=df_HU)

olog_PL <- polr(demand.ord, Hess = TRUE, method = "logistic", data=df_PL)

# Interaction models
demand.ord1 <- as.formula("demand_ord ~ aware_ct*exec_infl + df_support + perf_ct + copartisan + ideo10 + pol_interest + strong_lead + know_ct + age + gender")

olog_US1 <- polr(demand.ord1, Hess = TRUE, method = "logistic", data = df_US)

olog_DE1 <- polr(demand.ord1, Hess = TRUE, method = "logistic", data = df_DE)

olog_HU1 <- polr(demand.ord1, Hess = TRUE, method = "logistic", data = df_HU)

olog_PL1 <- polr(demand.ord1, Hess = TRUE, method = "logistic", data = df_PL)

## Sandwich doesn't work well with polr {MASS} objects
  # Function to compute HC3 standard errors by hand for polr objects
  # See https://library.virginia.edu/data/articles/understanding-robust-standard-errors
HC3_se <- function(polr) {
  library(sure)
  library(MASS)
  # Extract sigma from the polr model
  s2 <- sigma(polr)
  # Extract the model matrix
  X <- model.matrix(polr)
  # Number of observations
  n <- nrow(X)
  # Residuals from the polr model
  res <- sure::resids(polr)
  # Degrees of freedom
  df <- n - ncol(X) - 1 - 1
  # Hat matrix H
  H <- X %*% solve(t(X) %*% X) %*% t(X)
  # hat values
  hat_values <- diag(H)
  # HC1 calculation
  hc3 <- res^2/(1 - hat_values)^2
  # HC3 meat matrix
  hc3_meat <- hc3 * diag(n)
  # Variance-covariance matrix estimation
  vce_hc3 <- solve(t(X) %*% X) %*% (t(X) %*% hc3_meat %*% X) %*% solve(t(X) %*% X)
  # HC3 standard errors, excluding the intercept term
  hc3_se <- sqrt(diag(vce_hc3))[-1]
  return(hc3_se)
}

# HC3 standard errors using HC3_se function (see above)
  # additive models
se_US <- HC3_se(olog_US)
se_DE <- HC3_se(olog_DE)
se_HU <- HC3_se(olog_HU)
se_PL <- HC3_se(olog_PL)
  # interaction models
se_US1 <- HC3_se(olog_US1)
se_DE1 <- HC3_se(olog_DE1)
se_HU1 <- HC3_se(olog_HU1)
se_PL1 <- HC3_se(olog_PL1)

# Extract goodness-of-fit measures
res.var <- round(c(olog_US$deviance, olog_DE$deviance, olog_HU$deviance, olog_PL$deviance, 
                   olog_US1$deviance, olog_DE1$deviance, olog_HU1$deviance, olog_PL1$deviance),1)
aic <- round(c(AIC(olog_US), AIC(olog_DE), AIC(olog_HU), AIC(olog_PL), 
               AIC(olog_US1), AIC(olog_DE1), AIC(olog_HU1), AIC(olog_PL1)), 1)

# Table C3
labels <- c("Awareness", "Perceived Exec. Infl.", "Awareness $\\times$ Exec. Infl.", 
            "Diffuse Supp.", "Specific Support", "Gov. Supporter", "Ideology", 
            "Pol. Interest", "Strong Leader", "Knowledge Court", "Age", "Gender")

tabC3 <- stargazer(olog_US, olog_DE, olog_HU, olog_PL, 
          olog_US1, olog_DE1, olog_HU1, olog_PL1,
          se = list(se_US, se_DE, se_HU, se_PL,
                    se_US1, se_DE1, se_HU1, se_PL1),
          keep.stat = c("n","ll","aic","bic"),
          ord.intercepts = TRUE,
          star.char = c("+", "*", "**", "***"),
          star.cutoffs = c(.1, .05, .01, .001), 
          no.space=TRUE, 
          covariate.labels = labels,
          align=TRUE, table.placement = "H", 
          title="Determinants of Demand for Judicial Independence (Ordered Logit Models)", 
          dep.var.labels.include = FALSE, header = FALSE,
          dep.var.caption = "",
          column.labels = c("\\textbf{Additive Model}",
                            "\\textbf{Interaction Model}"),
          column.separate = c(4, 4),
          order = c(1:2, 12, 3:11,13), column.sep.width = "-32pt",
          add.lines = list(c("Res. Variance", res.var), 
                           c("AIC", aic)),
          notes.append=FALSE, notes.align = "l", notes.label = "",
          notes = c("\\footnotesize \\parbox[t]{\\textwidth}{\\textit{Note:} Robust (heteroskedasticity-consistent) standard errors in parentheses. \\hfill $^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001}"))
cat(tabC3, sep = '\n', file = "tab/JLC-tabC3.tex")

###########
# Variable statistics
# Appendix Tables B2-B5
###########

# Create tables to store variables statistics
sum_US <- data.frame("Variable" = c("Demand for Judicial Independence", "Awareness", 
                                    "Perceived Exec. Infl.", "Diffuse Supp.", 
                                    "Specific Support", "Gov. Supporter",
                                    "Ideology", "Pol. Interest", "Strong Leader", 
                                    "Knowledge Court", "Age", "Gender"), 
                     "Mean" = NA, "SD" = NA, "Min" = NA, "Max" = NA)
sum_DE <- sum_US
sum_HU <- sum_US
sum_PL <- sum_US

# vector of variables to summarize
vars <- c("demand", "aware_ct", "exec_infl", "df_support", "perf_ct", "copartisan", 
         "ideo10", "pol_interest", "strong_lead", "know_ct", "age", "gender")

# loop over country samples
for(i in seq_along(vars)){
  x <- vars[i]
  
  sum_US[i, 2] <- mean(df_US[[x]], na.rm = T)
  sum_DE[i, 2] <- mean(df_DE[[x]], na.rm = T)
  sum_HU[i, 2] <- mean(df_HU[[x]], na.rm = T)
  sum_PL[i, 2] <- mean(df_PL[[x]], na.rm = T)
  
  sum_US[i, 3] <- sd(df_US[[x]], na.rm = T)
  sum_DE[i, 3] <- sd(df_DE[[x]], na.rm = T)
  sum_HU[i, 3] <- sd(df_HU[[x]], na.rm = T)
  sum_PL[i, 3] <- sd(df_PL[[x]], na.rm = T)
  
  sum_US[i, 4] <- min(df_US[[x]], na.rm = T)
  sum_DE[i, 4] <- min(df_DE[[x]], na.rm = T)
  sum_HU[i, 4] <- min(df_HU[[x]], na.rm = T)
  sum_PL[i, 4] <- min(df_PL[[x]], na.rm = T)
  
  sum_US[i, 5] <- max(df_US[[x]], na.rm = T)
  sum_DE[i, 5] <- max(df_DE[[x]], na.rm = T)
  sum_HU[i, 5] <- max(df_HU[[x]], na.rm = T)
  sum_PL[i, 5] <- max(df_PL[[x]], na.rm = T)
}

# Convert to latex tables
table_US <- xtable(sum_US, caption = "Variables Statistics - US")
table_DE <- xtable(sum_DE, caption = "Variables Statistics - DE")
table_HU <- xtable(sum_HU, caption = "Variables Statistics - HU")
table_PL <- xtable(sum_PL, caption = "Variables Statistics - PL")

# save tables
sink("tab/JLC-tabB2.tex")
print(table_US, type = "latex", include.rownames = FALSE, 
      caption.placement = "top", table.placement = "H")
sink()

sink("tab/JLC-tabB3.tex")
print(table_DE, type = "latex", include.rownames = FALSE, 
      caption.placement = "top", table.placement = "H")
sink()

sink("tab/JLC-tabB4.tex")
print(table_HU, type = "latex", include.rownames = FALSE, 
      caption.placement = "top", table.placement = "H")
sink()

sink("tab/JLC-tabB5.tex")
print(table_PL, type = "latex", include.rownames = FALSE, 
      caption.placement = "top", table.placement = "H")
sink()

###########
# Diffuse Support Variable
# Appendix Table D1
###########

# Create table to store factor analysis statistics
tabD1 <- data.frame("col1" = c("Multiple R square of scores with factors", 
                        "Cronbach's $\\alpha$", "Factor 1 Eingenvalue", 
                        "Factor Loadings", " ", " "," ", " "),
                 "col2" = c(" ", " ", " ", 
                         "Do Away", "Jurisdiction Strip", 
                         "Mixed Up", "Remove Judges", "Less Independent"), 
                 "US" = NA, "DE" = NA, "HU" = NA, "PL" = NA)
tabD1

### Factor analysis

datasets <- list("US" = df_US, 
                 "DE" = df_DE, 
                 "HU" = df_HU, 
                 "PL" = df_PL)

for(i in seq_along(datasets)){
  country <- datasets[[i]]
  
  # Compute Cronbach's alpha
  alpha <- round(cronbachs_alpha(country[, c("doaway_ct", "strip_ct", "mixed_ct", 
                                             "judges_rm_ct", "less_ind_ct")]), 2)
  # Factor analysis
  fa <- fa(country[, c("doaway_ct", "strip_ct", "mixed_ct", 
                       "judges_rm_ct", "less_ind_ct")], 
     nfactors = 1, rotate = "none", fm="pa", scores = "regression")
  
  # Extract multiple R square of scores with factors
  R2 <- round(as.numeric(fa$R2),2)
  
  # Extract factor 1 eingenvalue
  ein1 <- round(fa[["Vaccounted"]][1,], 2)
  
  # Extract factor loadings
  loadings <- round(as.numeric(fa$loadings[,1]),2)
  
  # Store values
  tabD1[1,i+2] <- R2
  tabD1[2,i+2] <- alpha
  tabD1[3,i+2] <- ein1
  tabD1[4:8,i+2] <- loadings
  
  # Construct diffuse support scale using regression scoring
  scores <- fa$scores
  min <- min(scores, na.rm = TRUE)
  max <- max(scores, na.rm = TRUE)
  
  if(names(datasets)[i] == "US"){
    df_US$df_support2 <- as.numeric(
      (scores - min) / (max - min)
    )
  }
  if(names(datasets)[i] == "DE"){
    df_DE$df_support2 <- as.numeric(
      (scores - min) / (max - min)
    )
  }
  if(names(datasets)[i] == "HU"){
    df_HU$df_support2 <- as.numeric(
      (scores - min) / (max - min)
    )
  }
  if(names(datasets)[i] == "PL"){
    df_PL$df_support2 <- as.numeric(
      (scores - min) / (max - min)
    )
  }

}

tabD1

# Check replication of diffuse support scale is successful
cor(df_US$df_support, df_US$df_support2)
cor(df_DE$df_support, df_DE$df_support2)
cor(df_HU$df_support, df_HU$df_support2)
cor(df_PL$df_support, df_PL$df_support2)

### Appendix Table D1

# Convert to latex table
tabD1
table_D1 <- xtable(tabD1, caption = "Diffuse Support Scale Factor Analysis Solutions")
colnames(table_D1) <- c("", "", "US", "DE", "HU", "PL")
table_D1

note <- list()
note$pos <- list()
note$pos[[1]] <- c(nrow(table_D1))
note$command  <- c(paste("\\hline \n",
                            "\\multicolumn{6}{l}{\\footnotesize \\textit{Note:} Loadings from the first factor of an unrotated solution of a common factor analysis.}\n",
                            sep = ""))

# Save the complete output to a .tex file
sink("tab/JLC-tabD1.tex")
print(table_D1, type = "latex", include.rownames = FALSE, 
      sanitize.text.function = identity, 
      add.to.row = note,
      hline.after = c(-1, 0), 
      caption.placement = "top", table.placement = "H")
sink()


#################
# V-Dem Indicators
# Figures 1 and 2
#################

# Remove all objects
rm(list = ls())

vdem <- readRDS(file = "data/V-Dem-CY-Full+Others-v14.rds")

# high court independence (higher values = more independence)
vdem$vdem_ind <- vdem$v2juhcind
# reforms (higher values = less reforms undermining judicial power)
vdem$vdem_reform <- vdem$v2jureform*-(1)
# court packing (higher values = more packing)
vdem$vdem_pack <- vdem$v2jupack*(-1)
# gov attacks on judiciary (higher values = more attacks)
vdem$vdem_attack <- vdem$v2jupoatck*(-1)
# Liberal democracy index (higher values = more democratic)
vdem$vdem_libdem <- vdem$v2x_libdem


vdem_dat <- vdem[vdem$year %in% c(2003:2023) & 
                   vdem$country_name %in% c("United States of America", 
                                            "Germany", "Hungary", "Poland"),]
countries <- split(vdem_dat, vdem_dat$country_name)

pdf(file = "fig/JLC-fig1.pdf", width = 5, height = 3)
# .tiff file
#tiff(file = "fig/JLC-fig1.tif", width = 5, height = 3, units = "in", res = 300)
par(mar=c(2, 2, 0.5, 1), oma=c(0.5, 0.5, 0.5, 0.5))
plot(NULL, xlim = range(vdem_dat$year), ylim = c(0.3,0.9),
     xlab = "", ylab = "", main = "", 
     xaxt = "n", yaxt = "n")

axis(side = 1, at = c(seq(2003,2023,4)), labels = c(seq(2003,2023,4)))
axis(side = 2, at = c(seq(0.3,0.9, 0.1)), labels = FALSE)

text(x = 2001, # horizontal position of x labels
     y = seq(0.3,0.9, 0.1), 
     labels = seq(0.3,0.9, 0.1),
     adj = 0.5, # alignment (0=left, 0.5=center, 1=right)
     xpd = NA,
     srt = 90, # rotation of labels in degrees (ie, 45=45°)
     cex = 0.85)

lines(countries[['United States of America']]$year, 
      countries[['United States of America']]$vdem_libdem, 
      type = "l", col = "black", lwd = 2, lty = 1)

lines(countries[['Germany']]$year, countries[['Germany']]$vdem_libdem, 
      type = "l", col = "black", lwd = 2, lty = 2)

lines(countries[['Hungary']]$year, countries[['Hungary']]$vdem_libdem, 
      type = "l", col = "grey", lwd = 2, lty = 2)

lines(countries[['Poland']]$year, countries[['Poland']]$vdem_libdem, 
      type = "l", col = "grey", lwd = 2, lty = 1)

legend("bottomleft", legend = c("US", "Germany", "Hungary", "Poland"), 
       lty = c(1,2,2,1), lwd = 2, ncol = 2, bty = "n",
       cex = 1, col = c("black","black", "grey", "grey"))

dev.off()


pdf(file = "fig/JLC-fig2-court-pack.pdf", width = 5, height = 3)
# .tiff file
#tiff(file = "fig/JLC-fig2-court-pack.tif", width = 5, height = 3, units = "in", res = 300)
par(mar=c(2, 2, 2, 1), oma=c(0.5, 0.5, 0.5, 0.5))
plot(NULL, xlim = range(vdem_dat$year), ylim = c(min(vdem_dat$vdem_pack, na.rm = T), 
                                                 max(vdem_dat$vdem_pack, na.rm = T)+0.5),
     xlab = "", ylab = "", main = "Court Packing", 
     xaxt = "n")

axis(side = 1, at = c(seq(2003,2023,4)), labels = c(seq(2003,2023,4)))

lines(countries[['United States of America']]$year, 
      countries[['United States of America']]$vdem_pack, 
      type = "l", col = "black", lwd = 2, lty = 1)

lines(countries[['Germany']]$year, countries[['Germany']]$vdem_pack, 
      type = "l", col = "black", lwd = 2, lty = 2)

lines(countries[['Hungary']]$year, countries[['Hungary']]$vdem_pack, 
      type = "l", col = "grey", lwd = 2, lty = 2)

lines(countries[['Poland']]$year, countries[['Poland']]$vdem_pack, 
      type = "l", col = "grey", lwd = 2, lty = 1)

legend("topleft", legend = c("US", "Germany", "Hungary", "Poland"), 
       lty = c(1,2,2,1), lwd = 2, ncol = 2, bty = "n",
       cex = 1, col = c("black","black", "grey", "grey"))

dev.off()

pdf(file = "fig/JLC-fig2-attack.pdf", width = 5, height = 3)
# .tiff file
#tiff(file = "fig/JLC-fig2-attack.tif", width = 5, height = 3, units = "in", res = 300)
par(mar=c(2, 2, 2, 1), oma=c(0.5, 0.5, 0.5, 0.5))
plot(NULL, xlim = range(vdem_dat$year), ylim = c(min(vdem_dat$vdem_attack, na.rm = T), 
                                                 max(vdem_dat$vdem_attack, na.rm = T)+0.5),
     xlab = "", ylab = "", main = "Gov't. Attacks on Judiciary", 
     xaxt = "n")

axis(side = 1, at = c(seq(2003,2023,4)), labels = c(seq(2003,2023,4)))

lines(countries[['United States of America']]$year, 
      countries[['United States of America']]$vdem_attack, 
      type = "l", col = "black", lwd = 2, lty = 1)

lines(countries[['Germany']]$year, countries[['Germany']]$vdem_attack, 
      type = "l", col = "black", lwd = 2, lty = 2)

lines(countries[['Hungary']]$year, countries[['Hungary']]$vdem_attack, 
      type = "l", col = "grey", lwd = 2, lty = 2)

lines(countries[['Poland']]$year, countries[['Poland']]$vdem_attack, 
      type = "l", col = "grey", lwd = 2, lty = 1)

legend("topleft", legend = c("US", "Germany", "Hungary", "Poland"), 
       lty = c(1,2,2,1), lwd = 2, ncol = 2, bty = "n",
       cex = 1, col = c("black","black", "grey", "grey"))

dev.off()

pdf(file = "fig/JLC-fig2-jud-ind.pdf", width = 5, height = 3)
# .tiff file
#tiff(file = "fig/JLC-fig2-jud-ind.tif", width = 5, height = 3, units = "in", res = 300)
par(mar=c(2, 2, 2, 1), oma=c(0.5, 0.5, 0.5, 0.5))
plot(NULL, xlim = range(vdem_dat$year), ylim = c(min(vdem_dat$vdem_ind, na.rm = T), 
                                                 max(vdem_dat$vdem_ind, na.rm = T)+0.5),
     xlab = "", ylab = "", main = "High Court Independence", 
     xaxt = "n")

axis(side = 1, at = c(seq(2003,2023,4)), labels = c(seq(2003,2023,4)))

lines(countries[['United States of America']]$year, 
      countries[['United States of America']]$vdem_ind, 
      type = "l", col = "black", lwd = 2, lty = 1)

lines(countries[['Germany']]$year, countries[['Germany']]$vdem_ind, 
      type = "l", col = "black", lwd = 2, lty = 2)

lines(countries[['Hungary']]$year, countries[['Hungary']]$vdem_ind, 
      type = "l", col = "grey", lwd = 2, lty = 2)

lines(countries[['Poland']]$year, countries[['Poland']]$vdem_ind, 
      type = "l", col = "grey", lwd = 2, lty = 1)

legend("bottomleft", legend = c("US", "Germany", "Hungary", "Poland"), 
       lty = c(1,2,2,1), lwd = 2, ncol = 2, bty = "n",
       cex = 1, col = c("black","black", "grey", "grey"))

dev.off()

pdf(file = "fig/JLC-fig2-reform.pdf", width = 5, height = 3)
# .tiff file
#tiff(file = "fig/JLC-fig2-reform.tif", width = 5, height = 3, units = "in", res = 300)
par(mar=c(2, 2, 2, 1), oma=c(0.5, 0.5, 0.5, 0.5))
plot(NULL, xlim = range(vdem_dat$year), ylim = c(min(vdem_dat$vdem_reform, na.rm = T), 
                                                 max(vdem_dat$vdem_reform, na.rm = T)+0.5),
     xlab = "", ylab = "", main = "Court-Curbing Reforms", 
     xaxt = "n")

axis(side = 1, at = c(seq(2003,2023,4)), labels = c(seq(2003,2023,4)))

lines(countries[['United States of America']]$year, 
      countries[['United States of America']]$vdem_reform, 
      type = "l", col = "black", lwd = 2, lty = 1)

lines(countries[['Germany']]$year, countries[['Germany']]$vdem_reform, 
      type = "l", col = "black", lwd = 2, lty = 2)

lines(countries[['Hungary']]$year, countries[['Hungary']]$vdem_reform, 
      type = "l", col = "grey", lwd = 2, lty = 2)

lines(countries[['Poland']]$year, countries[['Poland']]$vdem_reform, 
      type = "l", col = "grey", lwd = 2, lty = 1)

legend("topleft", legend = c("US", "Germany", "Hungary", "Poland"), 
       lty = c(1,2,2,1), lwd = 2, ncol = 2, bty = "n",
       cex = 1, col = c("black","black", "grey", "grey"))

dev.off()
